
function [p q] = implicit(dt)
global U  h  g  p  q r  nt dx x h1 itermax;

dx1=1/dx;
dx2 = 1/(dx*dx);
dt2= 1/dt;
beta=0.5;
%---- Construct the spatial Discretization Matrix------%
A = zeros(x,x);
A(1,x) = -U*0.5*dx1;
A(1,2) =  U*0.5*dx1;
for k=2:x-1
    A(k,k-1) = -U*0.5*dx1;
    A(k,k+1) =  U*0.5*dx1;
end

A(x,x-1) = -U*0.5*(1/dx);
A(x,1)   =  U*0.5*(1/dx);

L=zeros(2*x,2*x);
L(x+1:2*x,x+1:2*x)=A;
L(1:x,1:x)=A;

L(1,2*x)   = h*dx2;
L(1,x+2)   = h*dx2;
L(x,x+1)   = h*dx2;
L(x,2*x-1) = h*dx2;
for k=1:x
    L(k,x+k)= -2*h*dx2;
    L(x+k,k)=g;
end

for k=2:x-1
    L(k,x+k-1) = h*dx2;
    L(k,x+k+1)= h*dx2;
end

b=h*h/3*dx2;
C=zeros(2*x,x);
C(1,x)= b;
C(1,2)= b;
C(1,1)= -2*b;
for k=2:x-1
    C(k,k-1)=b;
    C(k,k+1) =b;
    C (k,k) = -2*b;
end
C(x,x-1)= b;
C(x,1)= b;
C(x,x)= -2*b;

D = C(1:x,1:x);


n = -2/15*h*h*h*dx2;
E=zeros(x,x);
E(1,x)= n;
E(1,2)= n ;
E(1,1)= -2*n+h/3;

for k=2:x-1
    E(k,k-1)=n;
    E(k,k+1) =n;
    E (k,k) = -2*n+h/3;
end

E(x,x-1)= n;
E(x,1)= n;
E(x,x)= -2*n+h/3;


I =eye(2*x,2*x);
r= zeros(x,1);
r_new=r;
sol(1:x,1) =p;
sol(x+1:2*x,1) =q;
sol_new=sol;

for n=2:nt+1; 
    r=R_Solve(q);
%    rhs = (dt2*I - (1-beta)*L)*sol; 
    for iter =1:itermax        
        rhs = (dt2*I - (1-beta)*L)*sol + beta*C*(r_new+r);
        sol_new=(dt2*I+beta*L)\rhs;
        r_new = R_Solve(sol_new(x+1:2*x,1));
        
    end
    sol = sol_new;
    p = sol(1:x,1);
    q= sol(x+1:2*x,1);    
   
    if rem(n,100)==0
        refreshdata(h1,'caller') % Evaluate p in the function workspace
        drawnow
    end

end
display('Completed Successfully');